This Markdown encompass the code that was used to generate the gene expression related figures of the (manuscript)[https://doi.org/10.1101/2022.09.07.506982] titled: Single-cell transcriptomics of resected human traumatic brain injury tissues reveals acute activation of endogenous retroviruses in oligodendroglia.
set.seed(7)
library(ggplot2)
library(Seurat)
## Attaching SeuratObject
library(ggpubr)
library(biomaRt)
## Warning: package 'biomaRt' was built under R version 4.0.3
library(pheatmap)
library(RColorBrewer)
## Warning: package 'RColorBrewer' was built under R version 4.0.5
library(stringr)
library(openxlsx)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:biomaRt':
##
## select
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(clusterProfiler)
##
## clusterProfiler v4.5.0.992 For help: https://yulab-smu.top/biomedical-knowledge-mining-book/
##
## If you use clusterProfiler in published research, please cite:
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. The Innovation. 2021, 2(3):100141
##
## Attaching package: 'clusterProfiler'
## The following object is masked from 'package:biomaRt':
##
## select
## The following object is masked from 'package:stats':
##
## filter
library(org.Hs.eg.db)
## Loading required package: AnnotationDbi
## Warning: package 'AnnotationDbi' was built under R version 4.0.3
## Loading required package: stats4
## Loading required package: BiocGenerics
## Warning: package 'BiocGenerics' was built under R version 4.0.5
## Loading required package: parallel
##
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:parallel':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, parApply, parCapply, parLapply,
## parLapplyLB, parRapply, parSapply, parSapplyLB
## The following objects are masked from 'package:dplyr':
##
## combine, intersect, setdiff, union
## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, basename, cbind, colnames,
## dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
## grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
## order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
## rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
## union, unique, unsplit, which.max, which.min
## Loading required package: Biobase
## Warning: package 'Biobase' was built under R version 4.0.3
## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.
## Loading required package: IRanges
## Warning: package 'IRanges' was built under R version 4.0.3
## Loading required package: S4Vectors
## Warning: package 'S4Vectors' was built under R version 4.0.3
##
## Attaching package: 'S4Vectors'
## The following object is masked from 'package:clusterProfiler':
##
## rename
## The following objects are masked from 'package:dplyr':
##
## first, rename
## The following object is masked from 'package:base':
##
## expand.grid
##
## Attaching package: 'IRanges'
## The following object is masked from 'package:clusterProfiler':
##
## slice
## The following objects are masked from 'package:dplyr':
##
## collapse, desc, slice
##
## Attaching package: 'AnnotationDbi'
## The following object is masked from 'package:clusterProfiler':
##
## select
## The following object is masked from 'package:dplyr':
##
## select
##
library(AnnotationHub)
## Warning: package 'AnnotationHub' was built under R version 4.0.5
## Loading required package: BiocFileCache
## Warning: package 'BiocFileCache' was built under R version 4.0.3
## Loading required package: dbplyr
##
## Attaching package: 'dbplyr'
## The following objects are masked from 'package:dplyr':
##
## ident, sql
##
## Attaching package: 'AnnotationHub'
## The following object is masked from 'package:Biobase':
##
## cache
tbi <- readRDS("results-merged.rds")
tbi_seurat <- tbi[[1]]
tbi_res <- tbi[[2]]
getPalette = colorRampPalette(brewer.pal(9, "Greys"))
DimPlot(tbi_seurat, cols = getPalette(23)[5:20], label = T) + theme(legend.position = "None")
tbi_samples <- c("MJ_TBI_nr1_LT", "MJ_TBI_nr10_LT", "MJ_TBI_nr11_LFP", "MJ_TBI_nr16_RT", "MJ_TBI_nr19_RF", "MJ_TBI_nr2_LF", "MJ_TBI_nr20_RF", "MJ_TBI_nr21_RF", "MJ_TBI_nr3_RF", "MJ_TBI_nr6_RT", "MJ_TBI_nr7_LT", "MJ_TBI_nr8_RT")
control_samples <- c("TBI_HuBrainCTL_Nuclei501F_F", "TBI_HuBrainCTL_Nuclei501T_T", "TBI_HuBrainCTL_Nuclei502T_T", "TBI_HuBrainCTL_Nuclei529F_F", "TBI_HuBrainCTL_Nuclei529T_T")
DotPlot(tbi_seurat, features=c("MAP2", "DCX", "RBFOX1", "GRIN1", "HS3ST2", "GAD1", "GAD2", "CALB2", "CNR1", "GFAP", "AQP4", "GJA1", "SLC1A3", "PLP1", "MOG", "COL9A1", "VCAN", "PDGFRA", "P2RY12", "FYB1")) + coord_flip()+ theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + labs(x="")
FeaturePlot(tbi_seurat, c("RBFOX3", "GAD1", "GFAP", "PLP1", "VCAN", "FYB1"), ncol=3, cols = c("lightgrey", "red"))
celltypes <- c("5" = "Excitatory Neurons",
"0" = "Excitatory Neurons",
"6" = "Excitatory Neurons",
"2" = "Interneurons",
"4" = "Astrocytes",
"9" = "Interneurons",
"3" = "Interneurons",
"1" = "Oligodendrocytes",
"8" = "Microglia",
"7" = "OPC",
"14" = "Endothelial",
"13" = "Endothelial",
"10" = "?",
"11" = "?",
"12" = "?")
celltypes <- as.data.frame(celltypes)
celltypes$seurat_clusters <- rownames(celltypes)
clusters <- FetchData(tbi_seurat, "seurat_clusters")
clusters$barcode <- rownames(clusters)
clusters <- merge(clusters, celltypes, by="seurat_clusters")
colnames(clusters) <- c("seurat_clusters", "barcode", "celltypes")
rownames(clusters) <- clusters$barcode
tbi_seurat <- AddMetaData(tbi_seurat, metadata = clusters[,"celltypes", drop=F], col.name = "celltypes")
DimPlot(tbi_seurat, group.by = "celltypes",cols = c("Excitatory Neurons" = "#6d68a1",
"Interneurons" = "#312c6e",
"Oligodendrocytes" = "#579160",
"Astrocytes" = "#2b6b35",
"OPC" = "#104d1a",
"Microglia" = "#e8cc4d",
"Endothelial" = "#e3a94b",
"?" = "#c2831d"), label = T, label.box = T, label.color = "white") + ggtitle("")
clusters <- FetchData(tbi_seurat, "sample")
clusters$condition <- ifelse(clusters$sample %in% tbi_samples, "TBI", "Control")
tbi_seurat <- AddMetaData(tbi_seurat, metadata = clusters[,"condition", drop=F], col.name = "condition")
DimPlot(tbi_seurat, label = T, split.by = "condition")
sample_celltypes <- as.data.frame(table(tbi_seurat$sample, tbi_seurat$celltypes))
colnames(sample_celltypes) <- c("sample", "celltype", "Freq")
# Import metadata
tbi_metadata <- read.xlsx("samples_metadata.xlsx")
sample_celltypes <- merge(sample_celltypes, tbi_metadata, by.x="sample", by.y="ID")
sample_celltypes$key_condition <- paste(sample_celltypes$Condition, sample_celltypes$key, sep=" - ")
ggplot(sample_celltypes, aes(x=key_condition, y=Freq, fill=Condition)) + geom_bar(stat="identity", position = "dodge")+ facet_wrap(.~celltype, scales= "free", ncol=4) + theme_classic() +theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + labs(x="", y="Percentage of nuclei in cell type (%)")+ scale_fill_manual(values = c("#bab8b8", "#db8f8f"))
conditions_celltypes <- FetchData(tbi_seurat, c("condition", "celltypes"))
conditions_celltypes <- as.data.frame(table(conditions_celltypes$condition, conditions_celltypes$celltypes))
conditions_celltypes$Var2 <- factor(conditions_celltypes$Var2, levels = c("Excitatory Neurons", "Interneurons", "Oligodendrocytes", "Astrocytes", "OPC", "Microglia", "?", "Endothelial"))
total_num_cells_conditions <- table(FetchData(tbi_seurat, c("condition")))
conditions_celltypes$Freq <- conditions_celltypes$Freq/(ifelse(conditions_celltypes$Var1 == "Control", total_num_cells_conditions[["Control"]], total_num_cells_conditions[["TBI"]]))
ggplot(conditions_celltypes, aes(x=Var1, y=Freq, fill=Var2)) + geom_bar(stat="identity") +
theme_classic() + labs(x="", y="Ratio of cells", fill="") + ylim(c(0,1)) +
scale_fill_manual(values = c("Excitatory Neurons" = "#6d68a1",
"Interneurons" = "#312c6e",
"Oligodendrocytes" = "#579160",
"Astrocytes" = "#2b6b35",
"OPC" = "#104d1a",
"Microglia" = "#e8cc4d",
"Endothelial" = "#e3a94b",
"?" = "#c2831d")) + theme(text = element_text(size=15))
sample_celltypes <- FetchData(tbi_seurat, c("sample", "celltypes"))
sample_celltypes <- as.data.frame(table(sample_celltypes$sample, sample_celltypes$celltypes))
sample_celltypes$Var2 <- factor(sample_celltypes$Var2, levels = c("Excitatory Neurons", "Interneurons", "Oligodendrocytes", "Astrocytes", "OPC", "Microglia", "?", "Endothelial"))
ggplot(sample_celltypes, aes(x=Var1, y=Freq, fill=Var2)) + geom_bar(stat="identity", position = "fill") +
theme_classic() + labs(x="", y="Ratio of cells", fill="") +
scale_fill_manual(values = c("Excitatory Neurons" = "#6d68a1",
"Interneurons" = "#312c6e",
"Oligodendrocytes" = "#579160",
"Astrocytes" = "#2b6b35",
"OPC" = "#104d1a",
"Microglia" = "#e8cc4d",
"Endothelial" = "#e3a94b",
"?" = "#c2831d")) + theme(text = element_text(size=15)) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
ggplot(sample_celltypes, aes(x=Var1, y=Freq, fill=Var2)) + geom_bar(stat="identity") +
theme_classic() + labs(x="", y="Number of cells", fill="") +
scale_fill_manual(values = c("Excitatory Neurons" = "#6d68a1",
"Interneurons" = "#312c6e",
"Oligodendrocytes" = "#579160",
"Astrocytes" = "#2b6b35",
"OPC" = "#104d1a",
"Microglia" = "#e8cc4d",
"Endothelial" = "#e3a94b",
"?" = "#c2831d")) + theme(text = element_text(size=15)) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
s.genes <- cc.genes$s.genes
g2m.genes <- cc.genes$g2m.genes
tbi_seurat <- CellCycleScoring(tbi_seurat, s.features = s.genes, g2m.features = g2m.genes, set.ident = F, assay="RNA")
## Warning: The following features are not present in the object: MLF1IP, not
## searching for symbol synonyms
## Warning: The following features are not present in the object: FAM64A, HN1, not
## searching for symbol synonyms
tbi_seurat <- AddMetaData(tbi_seurat, metadata = ifelse(tbi_seurat$Phase == "G1", "non-cycling", "cycling"), col.name = "cellCycle")
cellcycle_scores <- FetchData(tbi_seurat, c("UMAP_1", "UMAP_2","condition", "G2M.Score", "S.Score"))
ggarrange(ggplot(cellcycle_scores, aes(x=UMAP_1, y=UMAP_2, colour=G2M.Score)) + geom_point(size=0.5) +
facet_wrap(.~condition) + theme_classic() + scale_colour_gradient(low="lightgrey", high = "red"),
ggplot(cellcycle_scores, aes(x=UMAP_1, y=UMAP_2, colour=S.Score)) + geom_point(size=0.5) +
facet_wrap(.~condition) + theme_classic() + scale_colour_gradient(low="lightgrey", high = "red"), ncol=1)
stat1_control <- FeaturePlot(subset(tbi_seurat, condition == "Control"), features = c("STAT1"))
stat1_tbi <- FeaturePlot(subset(tbi_seurat, condition == "TBI"), features = c("STAT1"))
stat1_control$data$condition <- "Healthy"
stat1_tbi$data$condition <- "TBI"
stat1 <- rbind(stat1_control$data,
stat1_tbi$data)
colnames(stat1)[4] <- "STAT1"
ggplot(stat1, aes(x=UMAP_1, y=UMAP_2, colour=(STAT1))) + geom_point(size=0.5, alpha=0.5) + facet_wrap(.~condition) + theme_classic() + scale_color_gradient(low = "lightgray", high = "red") + labs(colour="") + theme(text=element_text(size=15)) + ggtitle("STAT1")
stat2_control <- FeaturePlot(subset(tbi_seurat, condition == "Control"), features = c("STAT2"))
stat2_tbi <- FeaturePlot(subset(tbi_seurat, condition == "TBI"), features = c("STAT2"))
stat2_control$data$condition <- "Healthy"
stat2_tbi$data$condition <- "TBI"
stat2 <- rbind(stat2_control$data,
stat2_tbi$data)
colnames(stat2)[4] <- "STAT2"
ggplot(stat2, aes(x=UMAP_1, y=UMAP_2, colour=(STAT2))) + geom_point(size=0.5, alpha=0.5) + facet_wrap(.~condition) + theme_classic() + scale_color_gradient(low = "lightgray", high = "red") + labs(colour="") + theme(text=element_text(size=15)) + ggtitle("STAT2")
cd47_control <- FeaturePlot(subset(tbi_seurat, condition == "Control"), features = c("CD47"))
cd47_tbi <- FeaturePlot(subset(tbi_seurat, condition == "TBI"), features = c("CD47"))
cd47_control$data$condition <- "Healthy"
cd47_tbi$data$condition <- "TBI"
cd47 <- rbind(cd47_control$data,
cd47_tbi$data)
colnames(cd47)[4] <- "CD47"
ggplot(cd47, aes(x=UMAP_1, y=UMAP_2, colour=(CD47))) + geom_point(size=0.5, alpha=0.5) + facet_wrap(.~condition) + theme_classic() + scale_color_gradient(low = "lightgray", high = "red") + labs(colour="") + theme(text=element_text(size=15)) + ggtitle("CD47")
nlrc5_control <- FeaturePlot(subset(tbi_seurat, condition == "Control"), features = c("NLRC5"))
nlrc5_tbi <- FeaturePlot(subset(tbi_seurat, condition == "TBI"), features = c("NLRC5"))
nlrc5_control$data$condition <- "Healthy"
nlrc5_tbi$data$condition <- "TBI"
nlrc5 <- rbind(nlrc5_control$data,
nlrc5_tbi$data)
colnames(nlrc5)[4] <- "NLRC5"
ggplot(nlrc5, aes(x=UMAP_1, y=UMAP_2, colour=(NLRC5))) + geom_point(size=0.5, alpha=0.5) + facet_wrap(.~condition) + theme_classic() + scale_color_gradient(low = "lightgray", high = "red") + labs(colour="") + theme(text=element_text(size=15)) + ggtitle("NLRC5")
ciita_control <- FeaturePlot(subset(tbi_seurat, condition == "Control"), features = c("CIITA"))
ciita_tbi <- FeaturePlot(subset(tbi_seurat, condition == "TBI"), features = c("CIITA"))
ciita_control$data$condition <- "Healthy"
ciita_tbi$data$condition <- "TBI"
ciita <- rbind(ciita_control$data,
ciita_tbi$data)
colnames(ciita)[4] <- "CIITA"
ggplot(ciita, aes(x=UMAP_1, y=UMAP_2, colour=(CIITA))) + geom_point(size=0.5, alpha=0.5) + facet_wrap(.~condition) + theme_classic() + scale_color_gradient(low = "lightgray", high = "red") + labs(colour="") + theme(text=element_text(size=15)) + ggtitle("CIITA")
Cell type composition
rownames(tbi_metadata) <- tbi_metadata$ID
time <- tbi_metadata[,"Time_Post_Injury",drop=F]
time$sample <- rownames(time)
time_sample <- FetchData(tbi_seurat, "sample")
time_sample$barcode <- rownames(time_sample)
time_sample <- merge(time_sample, time, by="sample")
rownames(time_sample) <- time_sample$barcode
tbi_seurat <- AddMetaData(tbi_seurat, metadata = time_sample[,"Time_Post_Injury", drop=F], col.name = "time")
FeaturePlot(subset(tbi_seurat, condition == "TBI"), "time")
regions <- tbi_metadata[,"region",drop=F]
regions$sample <- rownames(regions)
regions_sample <- FetchData(tbi_seurat, "sample")
regions_sample$barcode <- rownames(regions_sample)
regions_sample <- merge(regions_sample, regions, by="sample")
rownames(regions_sample) <- regions_sample$barcode
tbi_seurat <- AddMetaData(tbi_seurat, metadata = regions_sample[,"region", drop=F], col.name = "region")
DimPlot(tbi_seurat, group.by = "region") + ggtitle("Regions")
set.seed(6)
tbi_seurat <- SetIdent(tbi_seurat, value="celltypes")
celltypes <- as.character(unique(Idents(tbi_seurat)))
deas_celltype <- sapply(as.character(celltypes), function(x) NULL)
for(celltype in as.character(celltypes)) deas_celltype[[celltype]] <- FindMarkers(tbi_seurat, group.by = "condition", ident.1 = "TBI", subset.ident = celltype)
for(celltype in as.character(celltypes)) deas_celltype[[celltype]]$gene <- rownames(deas_celltype[[celltype]])
for(celltype in as.character(celltypes)) rownames(deas_celltype[[celltype]]) <- NULL
for(celltype in as.character(celltypes)) deas_celltype[[celltype]] <- deas_celltype[[celltype]][which(deas_celltype[[celltype]]$p_val_adj < 0.01),]
gse_dotplot <- function(dea){
genelist <- bitr(rownames(dea), fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")
genelist <- merge(genelist, dea[,c("avg_log2FC"), drop=F], by.x="SYMBOL", by.y="row.names")
genelist <- genelist[order(genelist$avg_log2FC, decreasing = T),]
genelist_FC <- genelist$avg_log2FC
names(genelist_FC) <- genelist$ENTREZID
gse <- gseGO(geneList=genelist_FC,
ont ="ALL",
keyType = "ENTREZID",
minGSSize = 3,
maxGSSize = 800,
seed = T,
pvalueCutoff = 0.05,
verbose = TRUE,
OrgDb = org.Hs.eg.db,
pAdjustMethod = "BH")
return(gse)
}
gses_celltype <- sapply(as.character(celltypes), function(x) NULL)
for(celltype in as.character(celltypes)) rownames(deas_celltype[[celltype]]) <- deas_celltype[[celltype]]$gene
# Saved the results to xlsx at tables/GO_overrepresentation/seed_7/per_celltype/
# gses_celltype$`Excitatory Neurons` <- gse_dotplot(deas_celltype$`Excitatory Neurons`)
# gses_celltype$Interneurons <- gse_dotplot(deas_celltype$Interneurons)
# gses_celltype$Oligodendrocytes <- gse_dotplot(deas_celltype$Oligodendrocytes)
# gses_celltype$OPC <- gse_dotplot(dea=deas_celltype$OPC)
# gses_celltype$Microglia <- gse_dotplot(dea=deas_celltype$Microglia)
# gses_celltype$Astrocytes <- gse_dotplot(deas_celltype$Astrocytes)
# gses_celltype$Endothelial <- gse_dotplot(deas_celltype$Endothelial)
# gses_celltype$`?` <- gse_dotplot(deas_celltype$`?`)
Saved the results to xlsx at tables/GO_overrepresentation/seed_7/per_celltype/
gses_celltype_oligos <- read.xlsx("tables/GO_overrepresentation/seed_7/per_celltype/oligodendrocytes.xlsx")
# It really is just all the activated terms, but just to be clear
response_genes_oligos <- unique(unlist(str_split(gses_celltype_oligos[which(gses_celltype_oligos$Description %in% c("response to interferon-gamma", "cellular response to interferon-gamma", "cellular response to cytokine stimulus", "response to cytokine", "innate immune response", "interferon-gamma-mediated signaling pathway", "defense response")),"core_enrichment"], "/")))
response_genes_oligos <- bitr(response_genes_oligos, fromType="ENTREZID", toType="SYMBOL", OrgDb="org.Hs.eg.db", drop = F)
## 'select()' returned 1:1 mapping between keys and columns
tbi_seurat <- ScaleData(tbi_seurat, rownames(tbi_seurat))
## Centering and scaling data matrix
ensembl = useMart("ensembl",dataset="hsapiens_gene_ensembl")
gene.data <- getBM(attributes=c('hgnc_symbol'),
filters = 'go', values = gses_celltype_oligos$ID, mart = ensembl)
tmp <- DoHeatmap(tbi_seurat, features = gene.data$hgnc_symbol)
## Warning in DoHeatmap(tbi_seurat, features = gene.data$hgnc_symbol): The
## following features were omitted as they were not found in the scale.data slot
## for the RNA assay: EPRS1, CCL3L3, TRIM51G, IGHV3OR16-17, ILRUN, GARIN5A, TRIM75,
## RIGI, MATR3, TASL, STING1, C4B_2, SHFL, SEPTIN7, SEPTIN4, SEPTIN14, SEPTIN1,
## SEPTIN9, ATXN7, SEPTIN11, SEPTIN12, SEPTIN2, TUBB8B, SEPTIN3, SEPTIN6, SEPTIN8,
## SEPTIN10, SEPTIN5, BMERB1, DEFB109B
tmp_df <- reshape2::dcast(tmp$data[,1:3], formula= Feature~Cell, value.var = "Expression")
rownames(tmp_df) <- tmp_df$Feature
tmp_df <- tmp_df[,-1]
annotation_col <- FetchData(tbi_seurat, vars = c("condition", "celltypes"))
annotation_col <- annotation_col[order(annotation_col$celltypes, annotation_col$condition),,drop=F]
annotation_col <- annotation_col[which(rownames(annotation_col) %in% colnames(tmp_df)),, drop=F]
annotation_col$celltypes <- factor(annotation_col$celltypes, levels = c("Oligodendrocytes", "OPC", "Astrocytes", "Microglia", "Excitatory Neurons", "Interneurons", "Endothelial", "?"))
col.pal <- colorRampPalette(colors = c("lightgrey", "white", "#2d518a"))(50)
# Takes forever to plot...
# pheatmap(tmp_df[,rownames(annotation_col)], cluster_cols = F, show_rownames = F, cluster_rows = T, show_colnames = F, color = col.pal, annotation_col = annotation_col, gaps_col = c(which(!duplicated(paste(annotation_col$condition, annotation_col$celltypes)))[-1]-1, rep(which(!duplicated(annotation_col$celltypes))[-1]-1, 3)))
padj_genes <- rbind(deas_celltype$Oligodendrocytes["STAT1", "p_val_adj", drop=F],
deas_celltype$Oligodendrocytes["STAT2", "p_val_adj", drop=F],
deas_celltype$Oligodendrocytes["IRF1", "p_val_adj", drop=F],
deas_celltype$Oligodendrocytes["PARP14", "p_val_adj", drop=F],
deas_celltype$Oligodendrocytes["NLRC5", "p_val_adj", drop=F],
deas_celltype$Oligodendrocytes["CIITA", "p_val_adj", drop=F],
deas_celltype$Oligodendrocytes["IFI16", "p_val_adj", drop=F])
padj_genes$p_val_adj <- ifelse(padj_genes$p_val_adj > 0.05 | is.na(padj_genes$p_val_adj), "n.s.", format(padj_genes$p_val_adj, digits = 2))
library(ggpubr)
ggarrange(VlnPlot(subset(tbi_seurat, celltypes == "Oligodendrocytes"), features = c("STAT1"), pt.size = 0, group.by = "condition") + theme(legend.position = "None") + stat_summary(fun = mean, geom='point', size = 1, color="red") + geom_text(x=1.5, y=3, label=padj_genes["STAT1","p_val_adj"]),
VlnPlot(subset(tbi_seurat, celltypes == "Oligodendrocytes"), features = c("STAT2"), pt.size = 0, group.by = "condition") + theme(legend.position = "None") + stat_summary(fun = mean, geom='point', size = 1, color="red") + geom_text(x=1.5, y=3, label=padj_genes["STAT2","p_val_adj"]),
VlnPlot(subset(tbi_seurat, celltypes == "Oligodendrocytes"), features = c("IRF1"), pt.size = 0.3, group.by = "condition") + theme(legend.position = "None") + stat_summary(fun = mean, geom='point', size = 1, color="red") + geom_text(x=1.5, y=3, label=padj_genes["IRF1","p_val_adj"]),
VlnPlot(subset(tbi_seurat, celltypes == "Oligodendrocytes"), features = c("PARP14"), pt.size = 0, group.by = "condition") + theme(legend.position = "None") + stat_summary(fun = mean, geom='point', size = 1, color="red") + geom_text(x=1.5, y=3, label=padj_genes["PARP14","p_val_adj"]),
VlnPlot(subset(tbi_seurat, celltypes == "Oligodendrocytes"), features = c("NLRC5"), pt.size = 0, group.by = "condition") + theme(legend.position = "None") + stat_summary(fun = mean, geom='point', size = 1, color="red") + geom_text(x=1.5, y=3, label=padj_genes["NLRC5","p_val_adj"]),
VlnPlot(subset(tbi_seurat, celltypes == "Oligodendrocytes"), features = c("CIITA"), pt.size = 0.3, group.by = "condition") + theme(legend.position = "None") + stat_summary(fun = mean, geom='point', size = 1, color="red") + geom_text(x=1.5, y=3, label=padj_genes["CIITA","p_val_adj"]),
VlnPlot(subset(tbi_seurat, celltypes == "Oligodendrocytes"), features = c("IFI16"), pt.size = 0, group.by = "condition") + theme(legend.position = "None") + stat_summary(fun = mean, geom='point', size = 1, color="red")+ geom_text(x=1.5, y=3, label=padj_genes["IFI16","p_val_adj"]), ncol=4, nrow=2)
## Warning: Removed 2449 rows containing missing values (`geom_text()`).
tmp <- DoHeatmap(subset(tbi_seurat, celltypes == "Oligodendrocytes"), features = c("HLA-A", "HLA-B", "HLA-C", "HLA-E", "HLA-F", "HLA-DMA", "HLA-DOB", "HLA-DRA", "HLA-DQB1", "HLA-DPA1"))
tmp_df <- reshape2::dcast(tmp$data[,1:3], formula= Feature~Cell, value.var = "Expression")
rownames(tmp_df) <- tmp_df$Feature
tmp_df <- tmp_df[,-1]
annotation_col <- FetchData(tbi_seurat, vars = c("condition"))
annotation_col <- annotation_col[order(annotation_col$condition),,drop=F]
annotation_col <- annotation_col[which(rownames(annotation_col) %in% colnames(tmp_df)),, drop=F]
col.pal <- RColorBrewer::brewer.pal(9, "Reds")
pheatmap(tmp_df[c("HLA-A", "HLA-B", "HLA-C", "HLA-E", "HLA-F", "HLA-DMA", "HLA-DOB", "HLA-DRA", "HLA-DQB1", "HLA-DPA1"), rownames(annotation_col)], cluster_cols = F, cluster_rows = F, show_colnames = F, color = col.pal, annotation_col = annotation_col)
tmp <- DoHeatmap(subset(tbi_seurat, celltypes == "Oligodendrocytes"), features = c("MYO1E", "KIF5B", "PIP5K1A", "PSMB8", "PSMB9", "TAP1", "TAP2", "CD86", "CD80", "CD74", "CIITA", "NLRC5"))
tmp_df <- reshape2::dcast(tmp$data[,1:3], formula= Feature~Cell, value.var = "Expression")
rownames(tmp_df) <- tmp_df$Feature
tmp_df <- tmp_df[,-1]
annotation_col <- FetchData(tbi_seurat, vars = c("condition", "time"))
annotation_col <- annotation_col[order(annotation_col$condition, annotation_col$time),,drop=F]
annotation_col <- annotation_col[which(rownames(annotation_col) %in% colnames(tmp_df)),, drop=F]
annotation_col$time <- as.factor(annotation_col$time)
pheatmap(tmp_df[c("MYO1E", "KIF5B", "PIP5K1A", "PSMB8", "PSMB9", "TAP1", "TAP2", "CD86", "CD80", "CD74", "CIITA", "NLRC5"), rownames(annotation_col)], cluster_cols = F, cluster_rows = F, show_colnames = F, color = col.pal, annotation_col = annotation_col)
padj_genes <- rbind(deas_celltype$OPC["STAT1", "p_val_adj", drop=F],
deas_celltype$OPC["STAT2", "p_val_adj", drop=F],
deas_celltype$OPC["IRF3", "p_val_adj", drop=F],
deas_celltype$OPC["PARP14", "p_val_adj", drop=F],
deas_celltype$OPC["NLRC5", "p_val_adj", drop=F],
deas_celltype$OPC["CIITA", "p_val_adj", drop=F],
deas_celltype$OPC["IFI16", "p_val_adj", drop=F])
ggarrange(VlnPlot(subset(tbi_seurat, celltypes == "OPC"), features = c("STAT1"), pt.size = 0, group.by = "condition") + theme(legend.position = "None") + stat_summary(fun = mean, geom='point', size = 1, color="red") + geom_text(x=1.5, y=2, label=formatC(padj_genes["STAT1","p_val_adj"], format = "e", digits = 2)),
VlnPlot(subset(tbi_seurat, celltypes == "OPC"), features = c("STAT2"), pt.size = 0, group.by = "condition") + theme(legend.position = "None") + stat_summary(fun = mean, geom='point', size = 1, color="red") + geom_text(x=1.5, y=2, label=formatC(padj_genes["STAT2","p_val_adj"], format = "e", digits = 2)),
VlnPlot(subset(tbi_seurat, celltypes == "OPC"), features = c("IRF3"), pt.size = 0, group.by = "condition") + theme(legend.position = "None") + stat_summary(fun = mean, geom='point', size = 1, color="red") + geom_text(x=1.5, y=2, label=formatC(padj_genes["IRF3","p_val_adj"], format = "e", digits = 2)),
VlnPlot(subset(tbi_seurat, celltypes == "OPC"), features = c("PARP14"), pt.size = 0, group.by = "condition") + theme(legend.position = "None") + stat_summary(fun = mean, geom='point', size = 1, color="red") + geom_text(x=1.5, y=2, label=formatC(padj_genes["PARP14","p_val_adj"], format = "e", digits = 2)),
VlnPlot(subset(tbi_seurat, celltypes == "OPC"), features = c("NLRC5"), pt.size = 0, group.by = "condition") + theme(legend.position = "None") + stat_summary(fun = mean, geom='point', size = 1, color="red") + geom_text(x=1.5, y=2, label=formatC(padj_genes["NLRC5","p_val_adj"], format = "e", digits = 2)),
VlnPlot(subset(tbi_seurat, celltypes == "OPC"), features = c("CIITA"), pt.size = 0, group.by = "condition") + theme(legend.position = "None") + stat_summary(fun = mean, geom='point', size = 1, color="red") + geom_text(x=1.5, y=2, label=formatC(padj_genes["CIITA","p_val_adj"], format = "e", digits = 2)),
VlnPlot(subset(tbi_seurat, celltypes == "OPC"), features = c("IFI16"), pt.size = 0, group.by = "condition") + theme(legend.position = "None") + stat_summary(fun = mean, geom='point', size = 1, color="red")+ geom_text(x=1.5, y=2, label=formatC(padj_genes["IFI16","p_val_adj"], format = "e", digits = 2)), ncol=4, nrow=2)
tmp <- DoHeatmap(subset(tbi_seurat, celltypes == "OPC"), features = c("HLA-A", "HLA-B", "HLA-C", "HLA-E", "HLA-F", "HLA-DMA", "HLA-DOB", "HLA-DRA", "HLA-DQB1", "HLA-DPA1"))
tmp_df <- reshape2::dcast(tmp$data[,1:3], formula= Feature~Cell, value.var = "Expression")
rownames(tmp_df) <- tmp_df$Feature
tmp_df <- tmp_df[,-1]
annotation_col <- FetchData(tbi_seurat, vars = c("condition", "time"))
annotation_col <- annotation_col[order(annotation_col$condition, annotation_col$time),,drop=F]
annotation_col <- annotation_col[which(rownames(annotation_col) %in% colnames(tmp_df)),, drop=F]
annotation_col$time <- as.factor(annotation_col$time)
col.pal <- RColorBrewer::brewer.pal(9, "Reds")
pheatmap(tmp_df[c("HLA-A", "HLA-B", "HLA-C", "HLA-E", "HLA-F", "HLA-DMA", "HLA-DOB", "HLA-DRA", "HLA-DQB1", "HLA-DPA1"), rownames(annotation_col)], cluster_cols = F, cluster_rows = F, show_colnames = F, color = col.pal, annotation_col = annotation_col)
tmp <- DoHeatmap(subset(tbi_seurat, celltypes == "OPC"), features = c("MYO1E", "KIF5B", "PIP5K1A", "PSMB8", "PSMB9", "TAP1", "TAP2", "CD86", "CD80", "CD74", "CIITA", "NLRC5"))
tmp_df <- reshape2::dcast(tmp$data[,1:3], formula= Feature~Cell, value.var = "Expression")
rownames(tmp_df) <- tmp_df$Feature
tmp_df <- tmp_df[,-1]
annotation_col <- FetchData(tbi_seurat, vars = c("condition", "time"))
annotation_col <- annotation_col[order(annotation_col$condition, annotation_col$time),,drop=F]
annotation_col <- annotation_col[which(rownames(annotation_col) %in% colnames(tmp_df)),, drop=F]
annotation_col$time <- as.factor(annotation_col$time)
pheatmap(tmp_df[c("MYO1E", "KIF5B", "PIP5K1A", "PSMB8", "PSMB9", "TAP1", "TAP2", "CD86", "CD80", "CD74", "CIITA", "NLRC5"), rownames(annotation_col)], cluster_cols = F, cluster_rows = F, show_colnames = F, color = col.pal, annotation_col = annotation_col)